#include <iostream>
#include <stdio.h>
#include <math.h>
#include <vector>
#include <deque>
#include <sys/shm.h>

#include "Timer/timer.01.h"
#include "floor.h"

#define SQR(x) ((x)*(x))

typedef double T;
typedef unsigned long long int uint;


void *new_shm(int key, uint size, char* addr){
	shmdt(addr);
	int id = shmget(key, 1, 0666);
	shmctl(id, IPC_RMID, NULL); // and erase it
	id = shmget(key, size, IPC_CREAT | 0666); // set seg with 'key' and 'size'
	return shmat(id, NULL, 0);
}

T mU = 3.87e6; // kg
T xU = 415; // m
T tU = 1e-6; // sec

T g = -9.81*SQR(tU)/xU; 

uint N = 110;

T Cr = 0;

std::deque<T> xh; // history


int main() {

	FILE *fp = fopen("out.B.txt", "w");
	
	uint &n = *(uint*)new_shm(9110, sizeof(N), NULL); n = N;
	Floor *f = (Floor*)new_shm(9111, N*sizeof(Floor), NULL);
	
	Timer tr;
	tr.set();
	
	// initial conditions
	
	for (uint i=0; i<N; i++) {
		f[i].h = 1.0/(10*N);
		f[i].v = 0;
		f[i].a = g;
		f[i].m = 1; 
		f[i].fixed = true;
		f[i].color = Color(0.5,0.5,0.5,1);
	}
	
	for (uint i=0; i<N; i++)
		f[i].x = (i+1)/(long double)N *(1- 0.5*f[N-1].h);

	
	
	/*
	for (uint i=N-15; i<N; i++) {
		f[i].a = g;
		f[i].fixed = false;
		f[i].color = Color(1,0.5,0,1);
	}
	*/
	
	f[N-15].color = Color(1,0,0,1);
	f[N-1].fixed = false;
	f[N-1].m*=15;

	
	
	T Etot0(0);
	for (uint i=0; i<N; i++)
	Etot0 += .5*f[i].m*SQR(f[i].v) - f[i].m*f[i].x*g;
	

	T Eloss(0);
	
	uint n_col(0);
	
	T h(0); for (uint i=0; i<N; i++) h+=f[i].h; // total floors height

	// simulation	
	uint t=0;
	do {
						
		for (uint i=0; i<N; i++) {
			
			if (f[i].fixed) continue;
				
			// advance pos and vels
			f[i].x += f[i].v + .5*f[i].a;
			f[i].v += f[i].a;
			
			// check and treat collisions 
			for (uint j=0; j<i; j++) {
								
				T dx = f[i].x-f[j].x;
				T dx_min = 0.5*(f[i].h+f[j].h);
								
				if (fabs(dx)<dx_min) {
								
					f[i].x+=dx_min-dx;
					
					T dv = f[j].v-f[i].v;
					T fmfvi_plus_fmfvj = f[i].m*f[i].v + f[j].m*f[j].v;
					T fmi_plus_fmj = f[i].m + f[j].m;
					T vi = (Cr * f[j].m*+dv + fmfvi_plus_fmfvj) / fmi_plus_fmj ;
					T vj = (Cr * f[i].m*-dv + fmfvi_plus_fmfvj) / fmi_plus_fmj ;
					
					Eloss+= f[i].m*SQR(f[i].v)+f[j].m*SQR(f[j].v) - f[i].m*SQR(vi)-f[j].m*SQR(vj);
					
					f[i].v = vi;
					f[j].v = vj;
					
					f[j].fixed = false; 
					
					n_col++;				
				}
								
			}
			
			// ground condition 
			if (f[i].x-.5*f[i].h < 0) {
				Eloss+= f[i].m*SQR(f[i].v);
				f[i].x = .5*f[i].h;
				f[i].v*= 0;
				
			}							
			
		}
		
		// control
		if (t%uint(0.2/tU)==0) {
			
			printf("t=%5.2f s\t", t*tU );
			printf("v[N-1]=% 6.2f m/s | % 6.2f km/h \t", f[N-1].v*xU/tU, f[N-1].v*xU/tU/1000*3600 );
			T Etot(0); for (uint i=0; i<N; i++) Etot += 0.5*f[i].m*SQR(f[i].v) - f[i].m*f[i].x*g;	
			Etot+=0.5*Eloss;
			printf("dE = % .3e\n", (Etot-Etot0)/Etot0 );
			
		}
		
		uint odt = 0.2/tU; // sec
		if (t%odt==0) { 
			
			xh.push_back(f[N-1].x);
			
			if (xh.size()>3) {
				
				xh.pop_front();
				//double dt = 
				double v = (xh[2]-xh[0])/(2*odt);
				double a = (xh[2]-2*xh[1]+xh[0])/SQR(odt);
				
				fprintf(fp, "% 6.2f % 7.2f % 7.2f % 7.2f\n", (t-odt)*tU, xh[1]*xU, f[N-1].v*xU/tU, a*xU/SQR(tU));
			
			}
			
		}
		
		//if (t == uint(4.0/tU)) break;
		
		t++;
		// end condition
		if (f[N-1].x + .5*f[N-1].h - h  <= 0) break;
		
	} while (true);
					
	printf("t=%g [sec] | %u \n", t*tU, t );
	printf("Eloss=%g [J] \n", Eloss * mU*SQR(xU)/SQR(tU) );
	printf("n_col=%u \n", n_col );
	tr.get();
	tr.print();
	fclose(fp);
}

